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Closed Markovian networks of queues that have the product form 
in their stationary probability distributions are useful in the perform- 
ance evaluation and design of computer and telecommunication 
systems. Therefore, the efficient computation of the partition func- 
tion—the key element of the solution in product form— has attracted 
considerable effort. We present a new and broadly applicable method 
for calculating the partition function. This method can be applied to 
very large networks, which were previously computationally intrac- 
table. Most of the paper details applications of this approach to a 
network class which arose in modeling an interactive processor. We 
show that the partition function and derivatives such as mean values 
(response times, CPU utilizations, etc.) may be represented by integrals 
and their ratios. The integrands contain a parameter N which is 
large for large networks. Next, the classical techniques of asymptotic 
analysis are applied to derive three main power series expansions in 
descending powers ofN to correspond to normal, high, and very high 
usage. This work emphasizes multiple terms in the expansions for 
precision and error analyses. 

I. INTRODUCTION 

The theoretical results on the product form of the stationary distri- 
butions of large classes of Markovian queuing networks continue to 
have a profound influence on computer communications, computer 
systems analysis, and traffic theory. 1 " 4 These results make at least 
feasible the analysis and synthesis of the large systems of ever increas- 
ing complexity being considered in these areas. The subclass of closed 
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networks of queues is more difficult to analyze than the open networks 
because there is no stationary independence of the network nodes. 
However, the incentive for investigating the closed networks does exist 
since they have been used to model multiple-resource computer sys- 
tems, 2,5 multiprogrammed computer systems, 6-8 time-sharing, 2 and 
window flow control in computer communication networks; 910 net- 
works with external inputs subject to blocking require the analysis of 
a large number of closed networks. 1112 The closed network model that 
we shall use for illustrative purposes arose in the modeling of a central 
processor in a node of a computer network. This network is subject to 
a variety of processing demands. In recognition of the utility of closed 
networks, considerable research and commercial interest has been 
directed towards developing efficient procedures for computing the 
partition function (the normalizing constant), the only element of the 
product form solution requiring significant computation. 13-18 

However, as these existing recursive techniques are applied to the 
problems of particular interest in the Bell System, wherein the con- 
stituents of the closed chains are many and the number of chains are 
many, their shortcomings are observed to be severe in the amount of 
computing time and memory required and the accuracy attained. A 
more detailed account appears in Section 2.4 and Section IX. Briefly, 
the existing recursive techniques are largely ineffectual. 

We present a new way to view the problem, which surmounts many 
of the difficulties associated with large networks. The approach is 
broadly applicable — as indicated in Section X — even though the paper 
is a detailed account of applications to a specific class of closed 
networks. The new approach consists, first, of recognizing that the 
partition function may be written as an integral with a large parameter 
N present in the integrand to reflect the large size of the network. 
Next, the classical techniques of asymptotic analysis are applied to 
derive an asymptotic power series, typically in descending powers of 
N. The integrand will have fundamentally different properties in 
different ranges of the system parameters and this will require corre- 
spondingly different expansions. Thus, in this paper we develop three 
separate series expansions — Proposition 3 (Section IV), Proposition 12 
(Section VII), Proposition 17 (Section VIII) — each corresponding to a 
specific range of values of the usage parameter a. It is worth empha- 
sizing that, commensurate with an objective of providing solutions 
with any desired accuracy, we give procedures for generating multiple 
terms in the asymptotic expansions, not just the dominant term. In 
Section VIII, we unify the preceding results by giving a common 
expansion that holds uniformly in the system parameters. The uniform 
expansions introduce in a natural way the parabolic cylinder (or 
Weber) functions, a classical family of special functions with many 
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antecedents and ties with other special functions. 13 Besides duplicating 
the specialized expansions derived earlier, the uniform expansion 
makes available for use the many well-documented and tested expan- 
sions that are known in connection with parabolic cylinder functions. 

Section IX describes a user-oriented software package that has been 
written in C-language to implement the approach developed here. We 
supply results obtained by the package on four test problems that 
arose in analyzing performance of a Bell System project. Also reported 
are the results of a comparison with a well known, commercially 
marketed package that obtains solutions recursively. Our package is 
able to solve the large problems, which are well beyond the range of 
the other package, and, surprisingly, solve the small problems as well 
with errors that have small bounds. 

Section X provides the basis for extending the approach developed 
here to quite general multiprocessor, multidiscipline queuing networks. 
We show that for most networks that have been shown to have the 
product form in their stationary probability distribution, the partition 
function has an integral representation. The expansions appropriate 
for its computation are not considered here. 

Not surprisingly, the new representation of the partition function as 
an integral — the starting point of our computational procedures — may 
be exploited anew to derive analytical estimates and bounds of the 
quantities of interest, such as throughput, mean response time, etc. 
We demonstrate particularly in Section 5.3 that these formulas explic- 
itly exhibit the system parameters and as such are rather useful as 
design and synthesis aids. (The bounds are also useful as checks on 
the computational procedure.) Purely computational procedures by 
themselves do not yield this particular form of insight into system 
behavior. 

The asymptotic sequences used typically are power series in N~\ 
where recall N is the generic large parameter. 14 Thus, the number of 
terms required to achieve the desired accuracy decreases with increas- 
ing N. In contrast, with recursive solutions the computational com- 
plexity grows with the network size. Also, the asymptotic methods 
handle increased numbers of classes of constituents with little incre- 
mental difficulty, while the computational complexity in recursive 
methods grows geometrically. Thus, the contrasting techniques are 
not replacements for each other but complementary: loosely speaking, 
the recursions are most effective for smaller networks, while the 
asymptotic expansions are most effective for large networks. 

The contrasting behavior with respect to a large number of classes 
is of particular importance in computer communications where, as 
Reiser, 9 Schwartz, 3 and others have pointed out, traffic corresponding 
to each source-destination pair is treated as a separate class and the 
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network closure follows from the windows employed in flow control. 
Reiser has developed heuristics to cope with this situation. 

References 11, 15, and 16 also contain results pertaining to large 
networks. 

II. NETWORK MODEL AND KNOWN RESULTS 
2.1 Model 

In the model (see Fig. 1) each constituent, which may be thought of 
as a terminal or station, of the closed network spends alternating 
periods of time in the two nodes that constitute the network — the 
'think' node (also node 1) and the 4 cpu' node (node 2). The think time 
in each cycle for each constituent is an independent random variable 
with an exponential distribution. The time spent in the cpu node 
depends on many factors since it is here that there occurs interaction 
between constituents being serviced. We stipulate that the cpu disci- 
pline is 'processor sharing'* and that the desired service time (i.e., the 
time required to service the job if the entire cpu was dedicated to the 
job) is an independent exponentially distributed random variable. 1,4 
The 'think' node is thus an oo-server center and in the terminology of 
Ref. 4, nodes 1 and 2 are respectively Type 3 and 2 centers. 

We stipulate that there is sufficient statistical inhomogeneity 
amongst the constituents to justify the existence of several, say p, 
classes of constituents, with class i having Ki constituents, 1 < i < p. 
Statistical homogeneity applies within a class in that pa and pa will 
respectively denote the mean think time and the mean desired service 
time that are common to all in class i. The variations among these 
mean values may be quite substantial. 

Our involvement with this model arose while modeling behavior of 
traffic through a processing node of a computer network. The number 
of classes of constituents is at least five, namely, time sharing; inquiry/ 
response and data-base query; batch and remote job entry; messages 
and broadcast; data entry/collect and screen type jobs. The mean 
values {pij} are obtained from benchmark measurements. In another 
variation of this problem, a finer classification of constituents was 
considered. Our interest is in cases where the individual class popula- 
tions {Ki} extend to several hundred, while the number of classes is of 
the order of ten. 



* In the processor-sharing discipline there is no overt queuing because all, say n, jobs 
present in the node simultaneously receive service at 1/n times the rate given to a single 
job by the processor. Thus, the rate given to a single specific constituent fluctuates with 
time. This discipline is the limiting case of the round robin discipline as the time 
quantum given to each job becomes arbitrarily small. 
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Fig. l — (a) There are p classes of constituents — shown as terminals — with Kj constit- 
uents in class/ (b) Constituents spend alternate periods of time in the think node and 
processor-sharing cpu node. 



2.2 Product form solution 

If Nij is the stationary random variable denoting the number of 
constituents of class i in node j, and if m is the following stationary 
probability ir(nn, • • • , n p i\ «i2, • • • , I/.2) = Pr[N n = n u , • • • , N p i = 
rtpi\ Nja = /112, • • • , Np 2 = n p2 ], then it is known that with the left hand 
side abbreviated to «r(Hi, n2), 4 



7r(ni, n 2 ) = 



G(K) \.-i null \pi 



&S)l£Mjiins)- 



p p2 2 w 



a-1 n*2i 



(1) 



where G(K) = G(Ki, K 2 , • • • , K p ) is the normalization constant so 
chosen as to make the sum of all quantities in (1) equal to 1. Explicitly, 



G(K)= £ 

m„-0 



e in p " 






»>(aS> 



(2) 
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The function G( • ) defined on the integer lattice in R p is referred to as 
the partition function.* 

2.3 System performance 

A number of interrelated system perform ance me asures are obtained 
from the partition function. We start with iV,i(K), the mean number 
of constituents of class i in node 1 ('think'), and obtain directly 
from (2), 



JV.-i(K) = ptiG(K - e,)/G(K), (3) 

where e, is our notation for the vector with the i th component unity 
and all other components zero. Thus, G(K — e,) is the partition 
function associated with a new population with one less constituent in 
the I th class. From (3) and Little's theorem applied to class i and node 
1, we obtain for the throughput of constituents of class i, 

A,(K) = G(K- ei )/G(K). (4) 

The mean response time, i.e., time spent in the CPU in each cycle by 
class i constituents, is obtained again from an application of Little's 
theorem: 

ti(K) = KiG(K)/G(K - e.) - Pil . (5) 

Finally, the utilization of the cpu by constituents of class i, 

udK) A I . . . £ -^- „(n u n 2 ) - Pl2 G(K - e,)/G(K). (6) 
L n J2 

The important point to note is that all the mean values given in (3) 
through (6) are simply obtained from the knowledge of the partition 
function estimated at two neighboring points of the integer lattice in 
R p . As the above quantities are all closely related, we shall henceforth 
consider only the last, (u,(K)}. 

Higher moments of Nn, the random number of constituents of class 
i in node 1, may also be obtained from knowledge of the partition 
function: 



NUK) = {p?iG(K - 2e«) + Pil G(K - e,)}/G(K), (7) 

and, for i ¥>j, 

NnNMK) = Pi ipjiG(K - e, - e,)/G(K). (8) 

Of course, the moments of {.Afo} are easily derived from moments of 

{Nn}. 



* The distribution in (1) and (2) is also the stationary distribution of other networks. 
For example, if node 2 is first-come-first-served with class independent service rate 
1/P2, then (1) and (2) with pa ■ P2 is the solution. 
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2.4 Recursive solutions 

The above results explain why the problem of efficiently computing 
the partition function has excited so many researchers. 1517-21 For the 
problem at hand it is easy to arrive at the following recursion by 
established techniques: 

G(K) = i Py2 G(K - ej) + ft P 4\. (9) 

The boundary conditions are: G(K) = for K 2* 0. 

Observe that the partition functions themselves can be scaled on 
account of the linearity and the fact that only ratios of partition 
function values have physical content. By the same token, implemen- 
tations of (9), for large K, typically give rise to values which are either 
very small or very large leading to rather severe problems of overflow 
and underflow. Proper scaling is only marginally helpful. 

The main problems with implementing (9) are with respect to time, 
memory required during computation, and accuracy. A straightforward 
application of (9) would require an estimated K p iterations, where K 
is the generic class size. Similarly, the storage required would be 
approximately K*~\ Now these crude estimates can clearly be im- 
proved upon by simply pruning or avoiding the computation of inter- 
mediate lattice points, but this would be at the cost of increased 
algorithmic complexity, and the extent of the accrued benefits are not 
generally known. The underflow/overflow phenomenon that affects 
the accuracy of the scheme has already been commented on; a no less 
severe problem is accumulation of round-off errors in a large number 
of iterations. 

The recursive solution in (9) is one of several that can be generated 
by recently discovered techniques. However, all solutions that we are 
aware of are recursive and share to varying degrees the three broad 
categories of limitations just discussed. 

III. INTEGRAL REPRESENTATIONS 
3. 1 Partition function 

We start with Euler's integral 

n\= J e-'t n dt. (10) 

Jo 

Substituting for the middle term in braces in (2) we obtain* 



* If the range of integer subscripts is not stipulated explicitly, then the range is 
understood to be [1, p\. 
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J m p =o m,-o \ (Ay - m.j)\J \ mji ) 

= OlA^!}"J •"n(Wi + *A)^*fc (ID 

To obtain our final form for the partition function, let N be the 
generic large parameter to be associated with a large network. Also, 
for/ = 1, 2, • • • p, let 

Kj = ftjN, (12) 



rjA 



mean think time 
mean service time 



Pji 

Pj2 



class j 

= yjN. (13) 

The suggestion in this notation is that {/?,} and {yy} are 0(1), which is 
the situation to which our work is primarily directed. 

There is considerable latitude in selecting N. The following choice 
is certainly not essential but as it does ease some of the manipulations, 
we use it throughout the paper: 

N=(Y[rh 1/1Kj . (14) 

An implication of this choice of N is that 

E&logyy-O. (15) 

We return to (11) to observe 

G(K, "( n |)f e-'IIto + ')*'<** 

= { u ly) N ^ l j e ' Nz U(yj + ^) fijN dz. (16) 



Finally, we have after using (14), 
Proposition 1: 

G(K) = ^n^H e-™*dz, (17) 

where f(z) A z - I ft log(yy + z). D (18) 
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We shall see that typically there is no need to compute the term in 
braces in (17). 

3.2 Representations of mean values and some higher moments 

We have two options concerning representations of the physically 
interesting quantities in (3) through (7). Concerning ourselves with 
u,(K) given in (6), we may simply use the respective integral represen- 
tations for G(K - e.) and G(K) and obtain for i = 1, 2, • • • , p, 



ft e- NfU) dz 

MK>-£— ^s , (19) 

Y ' N I e- NfM dz 



where ft and /( • ) are defined analogously to N and f( • ) but for a 
network with one less constituent in the i th class. Obviously N and ft 
as well as /"( • ) and /( • ) are going to be close to each other, but the 
above option takes no further notice of this fact. 

In the contrasting option, we proceed as follows. Observe that 

1 G(K + e.) 
i*(K + *r l - TTTrn-^ (20) 

p,2 Cr(K) 

Now, from (16), 

G(K + e.) = jgfj (n |§) J to + ««"' II irj + «) *« 

where to obtain the last equation we have proceeded just as we did 
from (16) to (17). The point to note is that in the above expression N 
and /"(•) are identical to that used in the expression for G(K) in (17) 
and (18). Finally, from (20), for i = 1, 2, • • • , p 

Proposition 2: 

ze- Nf{z) dz 



^^- l -w^hT-^r lD <22) 



Notice that the ratio of integrals, the only quantity requiring significant 
computation, is common to all classes. 
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The higher moments given in (7) and (8) have similar representa- 
tions which are easy to derive. We give the form for one term that 
occurs in (7) which reveals the general pattern: 

1 G(K + 2e,-) 1 1 



.2 r<nr\ 7vr2 2 



pfi G(K) N< yi(p, + l/NHfii + 2/N) 



f ze- Nfl *dz [ z 2 e~ Nf{x) dz 

Jo , JO 



| e- Nm dz | e~ m *dz 

The feature to note is that the i th moment involves integrals of the 
form J? z J e~ Nf{z) dz, j - 0, 1, • • • , i. 

3.3 Properties of the function f 

We shall later need to recall certain properties of the function f(z) 
in (18), z > 0. As a consequence of (15), 

/(0) = 0. (24) 

Note 

/<«>( 2 )= 1 - £ ft/(y, + z) for *=1 

= (-l)'(i - 1)! 2 ft/(y, + 2)' for i > 1, (25) 

and the following alternating sign property that holds for i = 
2,3,..-, 

(-lYf U) (z) > 0, z > 0. (26) 

Also, for i = 2, 3, • • • , 

|/ W (z)|<=|/ W (0)|, 2^0. (27) 

As the second derivative is positive in [0, oo), the function is always 
convex. 

Let us view the function f(z) in the interval y m < z < oo, where 
y m A — min.y, (see Fig. 2). As the figure indicates, the derivative of the 
function tends to oo at both ends of the interval. Coupled with the 
convexity and the other already established facts, it follows that the 
function has a unique stationary point, a minimum, in (y m , oo). As in 
the figure, we let z denote this unique stationary point. 

This stationary point may be obtained as the unique real solution in 
(ym, °°) to the equation 

y-i Y> + 2 
Thus, the largest real root of a polynomial of order p gives i. 
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fU) 





Fig. 2— Sketches of the function f(z). 

It will be important for us to distinguish the cases z < and z > 0. 
The slope of the function f(z) at the origin effectively indicates which 
case holds: the first if the slope is positive and the second otherwise. 
Thus, a key parameter of the system is 

= I Kiln - l. (28) 



From the previous discussion we thus have that 

= f(0) if o<0 
= f(z) if a > 0. 



Min f(z) = f(0) if a<0 

22*) 



(29) 



Equation (29) summarizes the background on the stationary point 
as needed for most of the paper. For example, z is needed only if a > 
0. However, Section VIII is exceptional in that, while considering small 
possibly negative a, the corresponding stationary point z is required to 
be known. Note that as a — > 0, z - a/£ Pi/yh 

The parameter a, a > -1, is an indicator of the traffic intensity, 
with increasing a corresponding to higher traffic intensities. Our re- 
sults, theoretical and numerical, show that a > corresponds to heavy 
usage corresponding to close to 100 percent utilization of the CPU. 
'Normal' usage in large networks will certainly require a < and in all 
likelihood a will not be close to 0. For this reason, the most compre- 
hensive results given here are for the case a < 0. 

IV. ASYMPTOTIC EXPANSIONS FOR NORMAL TRAFFIC 
Throughout Section IV we shall consider a < 0. 
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4.1 Laplace' 8 method 

We shall first apply Laplace's method 14,22,23 to obtain asymptotic 
expansions (see Appendix A for notation) for the integral in the 
representation of the partition function in (17), and subsequently to 
the other integrals in (22) and (23). Laplace's method observes that 
when N is large, the minimum of Nf(z) at z = is very sharp and that 
the dominant contribution to the integral comes from the neighbor- 
hood of 0. 

In the integral 

/= J e - NfU) dz, (30) 

Jo 

let us change variables, 

u A f(z), (31) 

so that 

'-{"**(=)*■ (32 > 

To obtain dz/du, let us begin with the power series convergent in the 
neighborhood of 0, 

u = f(z) = I M (33) 

where fj = f {j) (0)/j\. Standard procedures allow the series to be 
reversed, i.e., give z as a power series in u. Appendix A elaborates on 
this procedure and explicitly gives the leading terms. Let the series 
obtained by this procedure be 



so that 



sAjOO-Zjcrf (34) 



^ = I aju J , (35) 

au j^o 



where a, — (j + l)gj+i. Substitution of this series in (32) yields 

/ ~ 2 a, I e~ Nu u j du 
i*» Jo 



and thus 
Proposition 3: 



I~?,T£k as AT^oo. □ (36) 
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The leading three terms of the series are obtained from 

ao - 1/fi, ai = -2/V/1; a 2 = 3(2/1- fxf»)/f\. (37) 

The proof that (36) is an asymptotic expansion is available from 
various sources. 14,22,23 Indeed an explicit proof is available from the 
bounds that we develop in the following section in the course of an 
error analysis. Nonetheless, we sketch a proof based on Watson's 
lemma 24 applied to the integral in (32) — the lemma is anyhow used 
later. Watson's lemma considers the integral Jo e~ Nu h(u)du in which 
(i) h(u) has a convergent power series expansion in the neighbor- 
hood of the origin, and 

(«) there exist constants c h c 2 such that | h(t) \ <cie C2t for t > 0, and 
asserts that an asymptotic expansion for the integral is obtained by 
replacing h(u) by its power series and integrating term by term. Thus, 
series reversion to obtain dz/du and a subsequent application of 
Watson's lemma gives the asymptotic expansion in (36). 

The reader will note for future reference that the asymptotic expan- 
sion in (36) may also be written as 

I~ Zg u) (0)/N j , (38) 

where g( • ) = f~\ • ), as follows from (33) and (34). 

Asymptotic expansions of integrals of the form I lk) = Jo z k e~ Nf(z) dz 
follow with only slight modifications. Thus, in lieu of (32) we now have 

/<*>= f e- Nu \z k ^\du. (39) 

In principle, it is straightforward to obtain a power series for (z k dz/ 
du), which is convergent in the neighborhood of from the power 
series in (34). Using the following as the defining relation* for the 
sequence {a}**}, 7 = 0, 1, 2, • • • , 

z" ? = ( 2 &A (lU+ Dgj+i» J ) = 1 «JV*. <40) 

du \j*l / \jso / /so 

the asymptotic expansion for the integral is obtained after term-by- 
term integration, giving 



Proposition 4: 



/(« = f z * e -m»dz ~ 2 (y + k)\ai k) /N j+k+x (41) 

Jo J-° 



as iV— >oo. D 



* In this notation, the sequence {<*>} in (35) and (36) is {af>). 
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Let us consider in greater detail the expansion of the integral J (1) = 
Jo ze~ Nf{z) dz, which will be needed if (22) is used to compute the mean 
values. We have derived the following recursive formula which effi- 
ciently generates {a}"} from the coefficients {aj 0) } that are needed 
anyhow for the expansion of /: for ally > 0, 

a? =i a?l k+x atjk. (42) 

In particular, the leading three terms of the expansion of I {X) are 
obtained from, 

ao 11 = l/ft a\ X) = -3/2//1; a? = 2(5/1 - 2/, /«//?. 

We have derived, but omit to give, a more general recursive for- 
mula — of which (42) is a special case — for generating {aj* +1) } from 

4.2 Asymptotic expansions for the utilizations 

As discussed so far, both of the two options stated in Section 3.2 for 
representing the mean values [see (19) and (22)] require the develop- 
ment of asymptotic expansions of two separate integrals and the 
subsequent computation of the ratio. Here we observe that a single 
asymptotic expansion exists for the mean values. The un- 
derlying reason is that the asymptotic sequence {N~ J }, j = 
0, 1, • • • , form a multiplicative asymptotic sequence. 14,22 

In particular, if as in (36) and (41), 

n (0) „<0) (0) 

P- Nflz) dz -£^- + £^- + 2' — +... 




and 



then 



1 



zp- NfM dz ~^-+ 2 + 2 
Ze dZ N* N* N 4 



ze~ NfU) dz 

e' NfM dz 








where the sequence {b,} is obtained by formal substitution. 

The following gives the leading terms for the utilizations derived by 
the above procedure. 
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Proposition 5: 

For i = 1, 2, • • • , p, as N -> «, 

{u,(K + e.)}" 1 ^ + 1/AO ~ y, + ^ + ^ + ^, (43) 

where 

A! = -l/a; A 2 = 4f 2 /a*; A 3 = -40/i/a 5 - 18/ 3 /a 4 . D 

In accordance with an earlier observation, all terms of the asymptotic 
series other than the first are independent of i. 

Proposition 5 contains a justification for treating the parameter a 
[see (28)] as an indicator of traffic intensity. Using only the dominant 
term, 

(uf(K + ef)} -1 ~r</A (44) 

and 

utilization of cpu = £ ut(K) ~ 1 + a. (45) 

A necessary caveat is that the above, as indeed all results in this 
section, has been derived for the assumption a < 0. 

Since the utilization as given by (45) can come close to unity even 
with a < 0, (45) justifies another earlier statement that for large 
networks normal usage will not extend beyond the range a < 0. 

V. ERROR ANALYSIS AND PERFORMANCE BOUNDS FOR NORMAL 
TRAFFIC 

Maintaining the restriction a < placed in the preceding section, we 
supplement in two directions the results obtained so far. First, we 
obtain essential results on the error incurred in truncating the expan- 
sions. These results containing information extending beyond what is 
required as proof of asymptotic expansions are needed for very prac- 
tical reasons, such as to know how many terms to use and, more 
importantly, to help define the regime of applicability. In the second 
part of the section, certain rather special properties of the functions in 
the representations are used to derive analytical bounds on the net- 
work performance measures. 

5. 1 Completely monotonic functions 

The following result on the function g( • ) = f~ l ( • ) is key to much of 
the error analysis: 

Proposition 6: 

(-iyg^(u) < for u > 0, j = 1, 2, 3, • • ■ D (46) 
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An inductive proof is given in Appendix B. By virtue of this result, 
g ll) ( • ) is a completely monotonic, or alternating, function (see Ref. 23). 
The importance of this property stems from the role of g a) (u) = dzj 
du in the integral representation (32). 

5.2 Error bounds 

In connection with the expansion of the integral / in Proposition 3, 
let R m denote the error that accrues if only the leading m terms are 
used, i.e., 

R m = I-"l jlaj 0) /N J+1 = I - J g°\0)/N J . (47) 

j-Q /-I 

Now by the mean value theorem, 25 for each u there is a £(m) in [0, u] 
such that 

m 

g M(u) = I g U) (0)u j - l /(j - 1)! + g (m+l) d)u m /m\ (48) 

y'-i 

On substitution in (32), 

m 

/= E g J (0)/N J + R m (49) 

7-1 



where 

1 



m! 



R m = —,\ e- Nu g {m+ »{£{u))u m du. (50) 



A simple corollary to (50) and Proposition 6 is 

Proposition 7: 

R m > if m is even, 
< if m is odd. □ 

This, of course, means that the terms in the asymptotic expansion 
alternate in sign and that the partial sums of the asymptotic expansion 
alternately over- and underestimate the true value of the integral in 
the following manner. 

Proposition 8: For m even, 



Jo 



m+l 

X g^W/N' < | e~ mz) dz < X £ w (0)/iV'. □ 

i-l 



The above is quite useful since in most situations the designer would 
much rather overestimate than underestimate a measure such as CPU 
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utilization. In this context, both the upper and lower bounds are 
required since ratios of integrals occur in the measures. 

An implication of Proposition 6 is that \g lm+1) U)\ < |tf (m+1) (0)|, 
£ > 0, which together with (50) gives 

Proposition 9: 

R m <g lm+1) (0)/N m+1 if mis even 
>£ (m+1) (0)/iV m+1 if mis odd. □ 

The above propositions thus state that the error is numerically less 
than the first neglected term of the series, and has the same sign. In 
particular, we have an explicit proof that (36) constitutes an asymptotic 
expansion. More generally, the above results show that, on account of 
the specialty of the integral, the main results that we require from an 
error analysis are already present in the expansion. 

It is useful to examine in detail g (4) (0) and thence the bound on R a : 

\R a \ < (6 4 |«| 2 + 10|a|&2& 3 + 156i)/(|a| 7 N 4 ), (51) 

where 6, - (l - 1)! 2 fij/yj. (A look at the proof in Appendix B will 
convince the reader of the presence of |a| 2m+1 N m+1 in the denominator 
of the bound for \R m \.) The bound does make the suggestion that in 
cases where a is so small [i.e., utilization is very large, see (45)] that 
a 2 N is itself small, then \R 3 \ is large. More generally, in the case of 
small a 2 N, the number of terms in the series requiring computation to 
meet specifications on the accuracy may be large. Later we return to 
consider this case further. 

5.3 Bounds on mean values 

The following bounds which supplement the computational proce- 
dures are presented to serve as aids in design and synthesis. For i = 1, 
2, •••,/>, 



Proposition 10: 



{ Ui (K + e,)}" 1 - 



Pi + 1/N 

1 + Vl + 8f 2 /a 2 N 
2|a|2V ' 



(Pi + 1/N) 

(52) 



a 



>i^i 1 , □ (53) 

2/ 2 V 1+ Vl + 16f 2 /(irct 2 N)J 
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Recall that 2/ 2 = / (2) (0) = £ Pj/yj while a (here a < 0) and N are as 
given in (28) and (14). 

To prove Proposition 10, we may use the representation in (22) for 
Ui in which case it remains to bound from above and below the pair of 
integrals appearing there. This is done in Appendix C by making use 
of the sign properties of the higher-order derivatives of the function 
/(■)• 

Notice that as N — » <», the upper bound in Proposition 10 on 
{k,-(K + e,)}" 1 approaches {y, - l/(aJV)}/(# + 1/2V), the sum of the 
leading two terms of the series for «, in (43). Also, as N — > », the 
expressions in (52) and (53) both approach 0. 

VI. EXPANSIONS FOR THE CASE a « O 

At the end of Section 5.2, we commented that the expansions given 
earlier may require a large number of terms in regions where a » 0. 
For this case, we here generate a somewhat different and more efficient 
series. We shall, however, be quite brief in our exposition because the 
uniform asymptotic expansions to be derived in Section VIII also allow 
appropriate expansions to be obtained. The ad hoc but direct treatment 
here is supplementary. 

We need notation that is specific to this section. Let 

Ki = biN + diJN 1 ■ , „ n /r.v 

n = a t (N + cVN) J I=1 ' 2 > •••>"• < 54 > 

Each of the variables {di} and c may be either positive or nega- 
tive. However, as our interest in this section is in a « and as a ~ 
2] hi/at — 1 as N — » oo, we require that 

Z bi/ai = 1. (55) 

In a computational procedure the above restriction poses no particular 
problem. 

Also, we shall mainly consider the integral 



JAdlT*) J e' 1 U (r< + t) K -dt. 
Jo 



(56) 



Comparison with (16) shows that the integral is related to the partition 
function thus: 

G(K) = (u^I. (57) 

As previously, the computation of the quantity in parentheses is not 
required to obtain the mean values. Our main result is 
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Proposition 11: As iV— » <», 

/- J Cy/V^"" 1 , (58) 

/-o 

where the sequence {c,} is given below. □ 

The proof of the proposition is in Appendix D. Here we comment 
on some features of the sequence: 

Cj = f e- A " 2/2 - B 'Hj(v)dp, (59) 

Jo 

where A = £ h/a 2 > and £ = c - £ di/«i and Hy(i») is a polynomial 
of degree 3/ in j* with coefficients that are fairly straightforward to 
obtain. The key point is that A, B as well as the coefficients of the 
polynomials are all 0(1), i.e., N does not enter into their definitions. 
Results given in Section VIII indicate how the coefficients c, given in 
(59) may be effectively computed. 
We give below the leading three polynomials in the sequence 

»M>: 

Ho(p) - 1, 

Hi(v) = c - (X di/a 2 )v 2 /2 + <Z bi/al)p 3 /3, 

H 2 (p) = -c(£ dt/aftv*/2 + (S di/al + c J bi/al)v>/3 

+ [(Zd i /a i ) 2 -2Zbi/ai]v 4 /8 

- (S <tya!)(Z bi/aftv 5 /6 + (J 6 f /af)V/18. (60) 

VII. ASYMPTOTIC EXPANSIONS FOR HEAVY TRAFFIC CONDITIONS 

Here we obtain asymptotic expansions for the basic integral in (17) 
and (18) for the case a > 0. For reasons similar to those discussed 
earlier for the case a < 0, the expansions to use for a ~ are in 
Sections VI and VIII. Hence, the expansions given below are for 
exceptionally heavy traffic conditions, where a is not only positive but 
also not close to 0. 

7.1 Laplace's method 

The key difference from the treatment in Section 4.1 is the presence 
of the singularity at z (see Fig. 2), which will be assumed to be known. 
For large N the dominant contribution to the integral 

1= \ e- Nfli) dz (61) 
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comes from the neighborhood of z. A Taylor series expansion around 
i gives 

f(z) - h - I fj(z ~ z)\ (62) 

where 

fj = f j) (z)/j\, 7 = 0,1,2,.... (63) 

In particular, fory > 2, 



£-<-D 7 



E #/(* + *)' 



>• (64) 



We make the following specific decomposition of J, which is conven- 
ient here and even more so in the error analysis to follow with 22 as in 
Fig. 2, 



/o/= f e - NlfU) -?o ] dz+ f 



+ f »-*«/«**«&{. (65) 

Jz., 

Consider the terms in turn starting with the middle term. If we let 

u A f(z) -? , 2 > 2, (66) 

and use the series in (62) for the right-hand side, then we may reverse 
the series, as discussed in Appendix A, to obtain 

2 - 2 A g(u) = l gju j/2 (67) 

7-1 

with the coefficient gj depending only on the coefficients f 2 , ft, • • • , 
fj. Now, 



J e -mm)-M dz= f ° e -f g m {u)du 



■i 



-fa 

e~ Nu X aju (j - l)/2 du, (68) 



where aj = (j + l)gj+i/2. The individual integrals in the sum will be 
recognized to be incomplete gamma functions. 

On returning to (65) and the first term in the right-hand side, we 
find by an identical argument that 
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f e -mm-k dz „ f e -Nu £ {-iya jU u - 1)/2 du. (69) 

Jo Jo >-° 

The two integrals in (68) and (69) may conveniently be combined to 
give 

H° °° r 

e m Im e -Nu j 2a v u J - l/2 du+ e- Nlfiz) - fo] dz. (70) 

Jo >-° Je 2 

At this stage, the following two approximations are made, with then- 
effects bounded in Section 7.2 in the course of the error analysis: the 
integration interval in the first term is extended to [0, <») and the 
second term is ignored. Nonetheless, the error analysis shows that 

I ~ e' Nh I e~ Nu I 2a 2j u J - 1/2 du, (71) 

Jo >"° 

giving Proposition 12. 
Proposition 12: As TV— > oo, 

I=[ e- m * ) dz~e- N ?'>i2rU+K)a2j/N j+1/2 . D (72) 
Jo J-° 

RecaU that r(V4) = V^ and for./' = 1, 2, • • • , 

Tij + %) - v£ n«-i a - %). 

We give the leading three coefficients: 

<*> = mitt*, 

a 2 = (Me)U5/i - 12/ 2 / 4 )//P, 

a 4 = (%56)(-64/i/ 6 + 224/I/3/5 + 112/1/2 

- 504/ 2 /i/ 4 + 28l/5)//P*. (73) 

The procedure for obtaining the asymptotic expansion for the 
integral 



Jo 



/<» = I ze ~ m *dz (74) 

is similar. Notice 



(7 (1> - zl) = e"** J (2 - i>" N( ' w "* } «fc; (75) 

Jo 
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we find it convenient to expand the integral on the right-hand side. 
The expression analogous to (71) is 

(7 (1) - zl) ~ e^ I e- Nu l 2a$- lU j - x/2 du t (76) 

Jo > =1 



where, for/ = 0, 1, 2, 



af -i's kgf-^igk, (77) 



which is to be compared with the expression for a, following (68). The 
following is a useful formula for efficiently generating the sequence 
{a 1 }*} from {a,}, which is needed anyhow for computing J: 

af - 2 l aj-kdk/ik + 1), j - 0, 1, 2, ... . (78) 

A-0 

Recognizing the gamma functions in (76) gives Proposition 13. 
Proposition 13: 

[W _ u m e -"?o j 2 r (j + i) a^/N i+l/2 (79) 

as N — » oo. D 

The asymptotic expansions in Propositions 12 and 13 may also be 
combined, as discussed earlier in Section 4.2 to yield an asymptotic 
expansion for the mean values. In particular, we obtain Proposi- 
tion 14. 

Proposition 14: As N -* oo, 

{u,(K + e.)}" 1 ^ + 1/N) ~ (y, + z)+^ + p 2 , 

where 

A, = -3/ 3 /4 ft, 

A 2 = (6 h h U ~ 15 fl h/S - 135 fl)/fl. □ (80) 

7.2 Error analysis 

The analysis to be presented supplements the result in Proposition 
12 and the error estimates to be given provide guidelines for the use of 
the expansions. The broad outline of the analysis have been suggested 
in Ref. 23. 

As in Section 5.2, let R n denote the error incurred when only n 
leading terms of the series in Proposition 12 is used, i.e., 
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R n = I - e~ N ' f ° V T (j + 1) 2a 2j /N J+1/2 . (81) 

For 7 we will use the expression given in (70) and decompose the error 
R n thus 

Rn = -€n.l(N) + €nAN) + € 3 (N), (82) 

where 

€n,i(N) = e~ Nfo [ e' Nu (V 2a 2j u J - l/2 \ du, (83) 

€„, 2 (N) = e - N/o I ° e~ Nu ( l 2o 2 >» / - 1/2 J d», (84) 

€3 (iV) = f e~ NM dz. (85) 



*2 



Thus, the three terms on the right-hand side of (82) respectively 
denote components arising from the extension of the integration inter- 
val from [0, -/o] to [0, oo) in (70) and (71), the use of only n leading 
terms from the infinite series in (71), and the neglect of the second 
term in the right-hand side of (70). Each component is now bounded. 
To bound e n ,i(N), we make use of known bounds on the incomplete 
gamma function: 23 

f" e- Nu u'- 1/2 du - T (j + \, - N? \/n' +1 ' 2 



e' 



"him L\\J +l ' 2 „ / . 1 



'W|/oD 



N J+i/2 lN\fo\-max(j-^,o\\ 



for n|/ | >max (>--,0J. 



Thus, 



\€n,i(N)\^ mf 2 . . 2 \a fl \\f Q \ j+y \ (86) 

iy\/o\ — On y=o 



where 8 n = max(n - 3/2, 0). 

In bounding e„, 2 (A0, we will postulate the existence of a finite valued 
On with the property that 

2 2a 2j u j ~ v2 < 1 2a 2n | u^'V-", < u < \ f \. (87) 

j-n 

This approach fails when o n is infinite but the characteristics of the 
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integral and the specific decomposition (65) that has been employed 
preclude this possibility. Let 

o n — max ^n(u), (88) 

M<l/bl 

where 

,7-1/2 

(89) 



\p n (u) = - In 
u 



2 202jU J 

j=n 



2a 2n u 



n- 1/2 



Small u is where (see Ref. 23) the danger of unbounded a„ is usually 
most manifest. However, as 



, . . 0.2n+2 f(l2n+4 0. 2/1+2 \ 

^ n (u) + I 2— )U+ ... 

0,2n \ &2n 0.2n J 



(90) 



when u —* + , no problem arises here. 

Using (87) in the defining expression for e„,2(A7) in (84) 

rl/ol 
| C „, 2 (A7) | < e- N f°\2 a2n \ e- (N -°» )u u n ~ 1/2 du 

Jo 

= e- N/o |2a 2/1 |r(n + %)/(N - a„)" +1/2 . (91) 

The last term to be considered from (82) is €3(N). We use the 
following property. 

f(z)>f'(z 2 )(z-z 2 ), z>2 2 , (92) 

which yields 

MW|£ Wte)' (93) 

a small quantity compared to the right-hand side of (91). This con- 
cludes the process of bounding the components of the error term R n . 
A corollary is the proof to Proposition 12. 

The bound in (91) is the largest component in the error bound. In 
examining (91), we observe that the condition in which the bound is 
large is when a 2n /N n is large. Now the expression for a 2n contains in 
the denominator a term/i n+1/2 , as (73) attests. Thus, when/iN is small, 
we expect the asymptotic expansions in Section 7.1 to be inefficient. 
We return to this case later in Section 8.2, where this as well as the 
similar difficulty encountered in Section 5.2 — where a was negative 
and small — is treated in a unified manner. 

VIII. UNIFORM EXPANSIONS 

This section has two objectives. The first is to show that there is a 
framework that unifies the expansions in Sections IV, VI, and VII. 
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This consists of showing that the integrals of interest may each be 
given by a common expansion valid uniformly for the entire range of 
values of the system parameters. These expansions turn out to be in 
parabolic cylinder (or Weber) functions. 13,25,26 The advantage derived 
is that these classical special functions have extensively documented 
expansions for the entire range of parameter values. 13,26 Indeed, using 
these expansions we sketch in Sections 8.3 and 8.4, at the cost of some 
duplication, derivations of the expansions obtained earlier in Sections 
IV and VII for a < and a > 0, respectively. The second objective is 
to derive a computationally efficient expansion for the case asO, i.e., 
where the stationary point z is very close to the boundary of the 
integration interval. The error analysis in Sections 5.2 and 7.2 has 
shown the need for a separate treatment. The expansion that is 
obtained for this case in Section 8.2 is obtained from an appropriate 
expansion of the parabolic cylinder functions. 

8. 1 Uniform expansions in parabolic cylinder functions 
Consider the integral 



-f 



e- Nf{z) dz (94) 



without restrictions on the parameter a. Following Friedman, 24 con- 
sider a change of variables from z to u given by 

v 2 - lav = f(z), (95) 

where a is a parameter of the transformation to be fixed later. The 
objective of the transformation is that the component of the integrand 
in braces below 

/= f e -JW«-** tej dv (96) 

satisfy the dual requirements of boundedness and a convergent power 
series, as required for an application of Watson's lemma (see Section 
4.1 following Proposition 3). [The reader may verify that the simpler 
transformation u = f(z) violates the boundedness requirement when- 
ever a > and z = z since f {l) (z) = 0.] For the transformation in (95), 

dz 2(v - a) 

Tv = ~T(z)- (97) 

This suggests the selection of the parameter a to be such that v = a 
when z = z, with the accompanying indeterminacy and the possibility 
of boundedness of dz/dv. This key clue does indeed give a unique map 
of the form in (95) with the desired properties, as summarized below. 

CLOSED MARKOVIAN NETWORKS 623 



Proposition 15: For z > 0, let 

v(z) = a + sgn(2 - z) Jf(z) - f(z), (98) 

where the constant a depends on all the system parameters: 



a=(sgna)V=7(i). (99) 

The transformation is monotonic, increasing and maps [0, <») to 
[0, oo). Also dv/dz is continuous and uniformly bounded. □ 

The transformation is used to derive for dz/dv a convergent power 
series £" C/r* in a neighborhood of the origin. This is achieved in three 
steps: 

(i) Use (98) to obtain 

v(z) = Aiz + A 2 z 2 + A 3 z 3 + • • • . (100) 

(ii) Reverse the series (see Appendix A) to obtain 

z(v) - Biv + B 2 v 2 + B 3 v 3 + • • • . (101) 

(Hi) Differentiate term by term to obtain 

-^ - Co + Cxv + C 2 v 2 + • • • , (102) 

dv 

where Cj = (j + l)B/+i. 

The reader may verify that the leading terms of the sequence (CJ) 
thus obtained are as follows [recall from (33) the definition fj = 

f U) (0)/M 

Co = 2a/a, 



a-* 



a \a J 



55V*-i 



(103) 



The above tacitly assumes that z and hence a have been evaluated. 

The power series expansion for dz/dv may now be substituted in 
(96) to yield 



J-° Jo 



e -N(v2-2av) v J d ^ (104) 



The integrals appearing above are simply related to the parabolic 
cylinder functions U(>, •); thus: for 7 = 0, 1, 2, • • • , 

I e-^-^v'dv = (m J w U[j + 1,- aJ2N). (105) 

Expansions for related integrals such as 
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■f 



/( i) = ze'Wdz (106) 



are only slightly different. Here the term dz/dv is replaced by zdz/dv, 
which has the power series expansion 



zdz/dv = £ Cfv j , (107) 



j-i 



where, see (101), 



Specifically, 



y-i 



Cf = I Bj- k C k . (108) 



«"-<¥)■ 



ci»-- — 



3 /2a \ /2a\ 



— ) fit - 1 



(109) 



a \<x J \oc J 

The following summarizes the expansions in parabolic cylinder 
functions of the two integrals of main interest, with the expansion 
valid for all a. 

Proposition 16: 

1= f e ~ Nf{z) dz= i Cj f e- N{v2 - 2av) v J dv 
Jo >=° Jo 

- e Na2/2 



hw$™ u ( j+1 *- a ' M )' <110) 

rd) = ( ze - N ^dz 

Jo 

-^ti^^H--^) (m) 

where the sequences (CJ) and {Cy'} are respectively as obtained by 
the procedures in (100) through (102) and (108). Specifically the 
leading terms are as given in (103) and (109). □ 

We should add that the above expansions are not strictly asymptotic 
expansions since the parabolic cylinder functions do not satisfy the 
requirements of asymptotic sequences for certain ranges of the param- 
eter a 2 N. u The interested reader will find in Ref. 27 a description of 
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the process for obtaining uniform asymptotic expansions of the inte- 
grals. However, we have not found it necessary to undertake the 
additional effort required to obtain the coefficients of the uniform 
asymptotic expansions. This is because for a 2 N small, the case treated 
below in Section 8.2 and of main interest, the functions U(j + V6, — 
a\/2N),j > 0, have all the desirable properties that are required of 
asymptotic sequences. 

A noteworthy property of the functions U( • , • ) that can be important 
in computations is that it satisfy the recursion 13 

xU(j + to, x) - U(j - % x) - (j + 1) U(j + % x). (112) 

8.2 Expansions for the case of a 2 N small 

To motivate the results to be given here, observe that the stationary 
point of the curve of f(z) (see Fig. 2) is close to when a is small 
(utilization of cpu is high), since 

z~a/(2f 2 ) as i->0. (113) 

On enquiring how the parameter a behaves for small z and a, we find 
from (99) that 

a =^T72 + °(« 2 ) ^ «-»0. (114) 

Thus, the case of small a, which we know from Section 5.2 requires 
special treatment, corresponds to small a. 

Small a is also implied by small fz and is therefore also of interest on 
account of the discussion in Section 7.2. This follows from 

a = i(/ 2 ) 1/2 + 0(£ 3/2 ) as i— 0. (115) 

For small a 2 N, it is known 13 that fory > 0, 

V7T 

where 

/J 2 ,/2 ~ y-1 

d* a Ta+7/2) u + 1)U + 3) ' ' ' u + ' " 1} ' * even 

-i iiu+i)/2> u+2)U+4) '" u+i - lh iodd - (117) 

In these relations, notice first in (116) the desirable presence of the 
powers of N in the denominator. Secondly, in connection with the 

626 THE BELL SYSTEM TECHNICAL JOURNAL, MAY-JUNE 1 981 



[pi j) + (aJN)n[ J) + (av / A0 2 J u2 y) +••■], (116) 



sequence {n\ j) }, observe that/^//i| y) = 2(y + i + l)/{(i + 1)(* + 2)}. 
Thus, for fixed y, the sequence converges rapidly to with increasing 
i. These observations state that when using the expression (116) in the 
expansions of Proposition 16, first, it is necessary to compute only 
up to small values of the index j and, secondly, with a>/N small, the 
computation of the bracketed quantity in (116) also needs very few 
terms. 
The following summarizes the important computational procedure 

described above. 

Proposition 17: For small a 2 N, 

where pi i] is in (117), and a, [Cj), {Cf} appear in Proposition 16. D 

8.3 Expansions for normal traffic 

For normal traffic, a < and consequently a < 0. If in addition, 
a <sc or, specifically, a 2 N »> 2 , then 13 



e -N(v*-2ao) v J dv 



fl -<w + <j*wr* (120) 



(-2aN) i+1 V Ma 2 N) 32(a 2 NY 
It can be shown after some manipulations, which we omit, that this 
expansion when substituted in Proposition 16 is identical to the main 
result of Section IV, namely, the expansion in Proposition 3. The 
bridging relation is 

(v-a) 2 = u + a 2 , (121) 

where u is the integration variable in Section IV [see (31)] and v is the 
similar variable in the uniform expansion [see (98)]. 

8.4 Expansions for heavy traffic 

For heavy traffic, a > 0, and therefore a > 0. When a 2 N » 1 as well, 
then 13 



I 



y/N \ 4a z N 



^-iKflK/-* ,... , im 
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Notice the departure from the expressions in (116) and (120) in the 
absence of N j in the denominator. 

It may again be established, although not in a simple manner, that 
the main result of Section VII, the expansion in Proposition 12, is also 
obtained by substituting the above expansion in Proposition 16. The 
bridging relation is 

(v - a) 2 = u, (123) 

where the variable u is as in (66). 

IX. COMPUTATIONAL NOTES 

The asymptotic expansions of integral representations of various 
quantities in the closed queueing networks discussed above, have been 
implemented as a user-oriented interactive package on Digital Equip- 
ment Corporation's VAX 11/780 operating under programmers work 
bench UNIX* system version 3. The package written in C-language 
has about 400 C-language statements and occupies about 60 Kbytes of 
storage. 

The number of classes and constituents! that the package can 
accommodate is so large that in effect no restriction is placed on these 
parameters. The other main features of the package are enumerated 
below. 

(i) The package is user oriented and easy to use. The user is 
prompted for relevant problem data. As this is supplied, validation and 
feasibility checks are made and the user informed of any errors. 

(ii) The output of the package includes all relevant statistics on 
each class (response time, utilization, etc.) including the percentage 
error incurred in the expansion. As an option to the last named, the 
user can display the various terms used in the expansion. 

(Hi) The package is partly adaptive in that it automatically detects 
the divergence of the asymptotic expansion and truncates the series at 
the point of divergence. 

(iv) Numerical stability is enforced by proper choice of N. 

Numerous computational experiments were performed to compare 
the efficacy of our package, called ASYM, to the current version of a 
popular commercially available package CADS. CADS is marketed by 
Information Research Associates and a version of it runs on a VAX 
11/780 operating under UNIX time sharing system. The test problems 
run on both these packages are real-world problems encountered in 
performance analysis of a Bell System project. 



* Trademark of Bell Laboratories. 

f In the computer science applications, the class sizes in closed networks (KJ) are 
referred to as degrees of multiprogramming. 
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The results of the experiments indicate that CADS is unable to 
solve our moderately sized, closed-queuing-network problems. Two 
features that accompany the breakdowns are noteworthy. One is 
numerical instability which manifests itself as overflows or underflows. 
Recent experience is that rescaling the service rates, an often-men- 
tioned device to combat the problem, does not help in substantially 
increasing the problem size. This accuracy problem is of course relieved 
with the use of more powerful floating-point machines like the CDC 
6600 (manufactured by Control Data Corporation) or IBM 3033. The 
other feature is the built-in limitations of the package. For instance, 
the current version of CADS limits the degree of multiprogramming 
for any one class to 100. There is also a limit of three classes. 

The above discussion is not intended to disparage the usefulness and 
success of CADS. CADS is extremely powerful in solving small 
queueing-network problems. Packages implementing recursions, as 
CADS does, and implementing expansions of integrals complement 
each other and, when integrated, will provide a powerful general- 
purpose package. 

The computational experiments (see Tables I to IV below) yielded 
the interesting fact that the asymptotic expansions are quite effective 
(i.e., yield a small percentage of errors) even for small problems, 
provided the right expansion is used. This in conjunction with the 
linear growth in computation time with increased number of classes 



Table I — Results for test problem 
No. 1* 



S?C? eof Utilization of CPU ^PUrSl? 
Multipro- Giyen b CADS CPU Given 

gramming J by AbYM 



10 

20 


0.0417 
0.083 


0.0414 
0.0829 


30 
40 


0.124 
0.166 


0.124 
0.165 


50 
60 


0.207 
0.249 


0.207 
0.248 


70 
80 


0.290 
Breakdown 


0.289 
0.331 


"90 
100 


Breakdown 
Breakdown 


0.372 
0.413 


110 
150 


Breakdown 
Breakdown 


0.618 
0.618 


200 


Breakdown 


0.819 



* Problem specification: 

No. of classes = 1 

Think time = 240 seconds 

Processing time = 1 second 
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Table IV(a) — Problem specification for 
test problem 4 



Class 


Service Rate 
of Infinite 
Server for 
This Class 


Service 

Rate of 

CPU for 

This Class 


Degree of 
Multipro- 
gramming 
for This 
Class 


1 

2 


0.0033 
0.033 


20.0 
2 


5 

5 


3 
4 


0.0033 
0.033 


4 
4 


6 

5 


5 
6 


0.033 
0.033 


4 
6 


5 

5 


7 
8 


0.033 
0.00033 

0.00055 
0.0033 


20 
0.6 

0.6 

0.6 


5 
5 


9 
10 


5 

5 


11 
12 


0.00033 
0.00055 


0.2 
0.2 


5 
5 


13 

14 


0.0003 
0.033 


0.2 
1 


5 
5 


15 
16 


0.00033 
0.00055 


1 
1 


5 
5 


17 


0.0003 


1 


5 



Table IV(b) — Results of test problem 4 
output by ASYM 



Class 


Response 

Time 
(seconds) 


Utilization 


% Error in 
Utilization 


1 

2 


0.09 
0.834 


0.008 
0.080 


0.05 
0.045 


3 

4 


0.43 
0.42 


0.004 
0.046 


0.05 
0.048 


6 

6 


0.42 
0.28 


0.040 
0.027 


0.049 
0.049 


7 
8 


0.09 
2.85 


0.008 
0.003 


0.05 
0.06 


9 
10 


2.9 
2.82 


0.005 
0.027 


0.05 
0.05 


11 
12 


8.530 
8.523 


0.008 
0.013 


0.05 
0.05 


13 
14 


8.540 
1.63 


0.007 
0.156 


0.05 
0.04 


15 

16 


1.71 
1.72 


0.002 
0.003 


0.05 
0.05 


17 


1.71 


0.001 


0.05 
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makes the use of the asymptotic expansions attractive even in cases 
where the recursive implementation does not break down. 

Tables I and II display the results of problems run both on CADS 
and ASYM. Tables III and IV show the results on two large problems 
that were not admitted by CADS, but were solved with good accuracy 
by ASYM. 

It may be observed from Table I that the results from ASYM and 
CADS agree rather well, even for small degrees of multiprogramming. 
In the cases solved by both CADS and ASYM, a was small and N 
(about 240) large. Observe also that CADS was unable to solve cases 
with Ki larger than 70, even though higher values of K\ correspond to 
quite low usage of the cpu and are quite interesting. 

Table II shows that CADS also broke down on a relatively small 
problem with 2 classes and 60 customers in each class. 

Tables III and IV display the output given by ASYM for the large 
problems. Table III corresponds to a problem where the total degree 
of multiprogramming is 3500. For Table IV the total degree of multi- 
programming is 85, but there are 17 classes. As shown, the percentage 
error in both cases is well within an acceptable range. 

In conclusion, the computational experiments suggest that the ap- 
proach based on expansions of integral representations is robust, 
computationally fast, and to be recommended for a variety of problems, 
large and small. 

X. GENERALIZATIONS: INTEGRALS IN NETWORKS WITH MANY 
PROCESSORS 

This section shows that for a large class of closed Markovian 
networks, the partition functions possess simple integral representa- 
tions. This result provides the basis for future work on the computa- 
tions of the integrals from expansions rather like those given in this 
paper. The above-mentioned class of networks allows an arbitrary 
number of service centers with flexibility in operating disciplines. It is 
in fact the same class of networks shown by Baskett et al., in Ref. 4, 
that has the product form in the stationary distributions, except that 
we do not allow the service rate in Type 1 centers to depend on the 
number in queue. (To some extent this is only for convenience because 
for some specific and interesting dependencies, we have obtained the 
integral representations.) 

The representation of the partition function that is obtained is as a 
multiple integral, i.e., as an integral in Euclidean space of dimension q 
where q is the number of queuing centers, which are the centers of 
Type 1, 2, and 4 in the notation of Ref. 4. However, in spite of the 
complexity of the partition function, the form of the integrand is 
remarkable for its simplicity. 

632 THE BELL SYSTEM TECHNICAL JOURNAL, MAY-JUNE 1 981 



10.1 Background on product form solutions 

As previously, we let the number of classes of constituents be p. We 
henceforth consistently index the classes by the symbol y ; when the 
index for summation or multiplication is omitted, it is understood that 
the missing index is y , where 1 <y < p. A total of s service centers are 
allowed. We will find it natural to distinguish the centers of Types 1, 
2, and 4, which have queuing, from the remaining centers of Type 3, 
which do not. Thus, centers 1 through q will be the queuing centers 
while (q + 1) through s will be the Type 3 centers, which have also 
been called think nodes and infinite server nodes. 

Let the 

stationary state probability = 7r(y x , y 2 , • • • , y 8 ), (124) 

y. 4 (ni„ n 2 „ • • • , n pi ), 1 < i < s, 
riji A number of class— y jobs in center i. 

The well-known results on Markovian closed queuing networks with 
product form solutions may be given in the following form: 1,4 

7r(yi, ...,y.)-ill "to). < 125 > 

tri-i 

where *(y<) = Q>„)! II fej > 1 ^ * =S q, 

-nf^j). iq + D*i*8. (126) 

In the above formulas we have taken into account the previously 
stated assumption, namely, for the first-come-first-served discipline in 
Type 1 centers the service rate is independent of the number of jobs 
in queue. Also, in (126), 

expected number of visits of classy jobs to center i 
service rate of classy jobs in center i 

where the numerator is obtained from the given routing matrix by 
solving for the eigenvector corresponding to the eigenvalue at 1. 
In (125) G is, of course, the partition function and it is explicitly 

G= I ••• I n^y,), (127) 

l'n,-/r, l'n p =Ap «-l 

where we have written l'n, for £*-i n,, and the condition 1'nj = Kj to 
indicate the conservation of jobs in each class. Thus, 
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G = l 






n in 

i—q+l 



rijil 



(128) 



10.2 Integral representation 

Using Euler's integral, see (10), 

- f ■ I exp( - l uA £ 



1 

l'n„-K„ 



n(n» 



i-i 



1=9+1 



n in?; 



nn\ 



dui">du q . (129) 



Now by the multinomial theorem, 



•I! X P/«"« + S P.. dm • • • dug. (130) 

V_l '-7+1 / 

It is noteworthy, but not surprising, that the parameters p# for all the 
Type 3 centers appear lumped together. 

We now introduce the large parameter N and define 



h = % lsy*A 



exactly as in (12). However, we define 
P/i 



Yfi = 



N, 



1 <></>, l<i<g, 



(131) 



(132) 



£ P/m 

m— 9+1 



which is reciprocal to the natural extension of the parameters {yj) 
defined in (13). In common with Section 3.1, the suggestion in the 
notation is that in the generic large network {y yj } and {/?,} are O(l). A 
tacit assumption being made is that all job classes are routed through 
at least one Type 3 (infinite server) center. 

On substituting (131) and (132) in (130) and after a change of 
variables we obtain a form for the partition function, which is sum- 
marized below. 

Proposition 18: 



g = ( n q n 

y-i 



i pji) m 

i-q+l 



e~ mz) dz, (133) 
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where 

Z - [Zi, 2 2 , • • • , Zq]', 
/■(Z) = 1'Z- i & log(l + r/z), 

1 = [1, 1, ... , 1]', 

^ = [7,1, Y>2, • • • , Yy 9 ]', l^j^P- □ 

As before, the term in brackets in (133) will typically not be required 
to be computed. 

Future work will consider the expansions appropriate for the com- 
putation of the integral in (133). 

APPENDIX A 

Notation for Asymptotic Expansion, Series Reversion 

Asymptotic expansion: A series Y)=o a j/ NJ *■ ^d to De an asymp- 
totic expansion 14,22,23 of a function F(N) if 

F(N) - X aj/W = 0(N~ n ) as N -► oo 
>-o 

for every n = 1, 2, • • • . We write 

F(N) ~ £ o,/M 

7-0 

The series itself may be either convergent or divergent. 

Series reversion: Ifu = f(z),iio = f(z ), f'(z ) ^ 0, then by Lagrange's 
expansion 13,25 



Z = Zo + 2* M 



d j ~ l I z-zo 



dz J - l \f(z)-Uo 



(134) 



In particular, if /"(•) is specified in a Taylor series, the above expansion 
identifies the coefficients { gj) in the reversed power series z - z = 
27_i gj\u - Uo) J . We specifically identify the leading coefficients of two 
reversed series that have been used in Section 4.1 and Section 7.1. 
If 

u = fxz + f 2 z 2 + • • • , (135) 

then z = giu + g 2 u 2 + ■ • • , (136) 

where 

figi = 1, 
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n& = -h 

f\g3 = 2f 2 -hfz, 

flg< = 5f 1 f 2 f3-tff4-5fl, 

figs = Sfif 2 f4 + 3/f/l + 14/1 - fxh ~ 21/i/i/a. 
Similarly, if 

u = f 2 z 2 + fe 3 + f 4 z 4 + ... , (137) 

then z = giu 1/2 + g 2 u + gau 3/2 + • • • , (138) 

where 

tPm = i, 

fig* - -/a/2, 
/P* = (5/1 - 4/ 2 / 4 )/8, 
/!*- (8/2/3/4 -/U-2/!)/2, 

/2 3/2 ^6 - (-64/i/e + 224/i/ 3 / 5 + 112/f/l - 504/ 2 /i/ 4 + 231/3)/128. 

APPENDIX B 

Proof of Proposition 6 

Before giving the proof, it is worthwhile to generate expressions for 
the leading derivatives of g(.). In the notation of Section 4.1, u = f(z) 
and 2 =#(m), 

^ 2) ( M ) = -/ (2) (2)/(/ (1) (2)) 3 , 

S (3, (a) = [-/ (3) U)/ (1, (2) + 3(/ (2, (z)) 2 ]/(/ (1) (z)) 5 , 

5 <4) (") = [-/ (4) U)(/ (1, (2)) 2 + 10/ (1) (*)/ (2, (2)/ (3) ( 2 ) 

-15(/ (2) (z)) 3 ]/(/ (1) (2)) 7 . (139) 

For notational convenience, let y, and <£, denote, in this appendix only, 
g li) (u) and f {i) (z) respectively. Recall that <fn > and that (-l)ty, > 
for i = 2, 3, • • • . We will show by induction that (— l) n y„ < 0, n = 
1,2,3,.... 
Let the induction hypothesis be the following 

n even: tf- x y H (*i)"(X) + foi) 2 * 1 ^) • • • , (140) 

where < i; max [2i, 2i + 1] < 2n — 1; X and Y are arbitrary products 
of <J>2, <f>3, • • • , </>n with arbitrary positive numerical weight and the 
property that X > 0, Y < for all z > 0. 
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n odd: tf-V = • • • + (<tn) 2i (U) - (<f>i)* i+1 (V) • • • , (141) 
where U and V are like X and Y including that U > 0, V < for all 

2>0. 

For the proof, take the case of n even, first. A key observation 
is that since 91 does not occur in either X or Y, dX/dz < and 
dY/dz > 0. That is, in common with the functions 92, 93, ••- the 
derivative has opposite sign from the function. Also, from (140) 

y n = ■ • • — 2n-2i-l """ .2/1-21-2 " ' i4 ^' 

9l 9l 

Differentiate with respect to u, 

1 (-X(2n - 2i - D92 X' 



Yn+l — • • • I ,2n-2« r .2n-2i-l 

9l \ <Pl <P1 

1 / -Y(2/i - 2i - 2)<h Y' 

"*" T" I ^2n-2t-l + . 2n-2i-2 ' ' 

<f>l \ 91 91 t 

Hence, 

9 2 " + Vi + (fa) 2 '(X(2n - 2i - D92) 

- (pi) 2,+1 (X' + Y(2n - 2i - 2)92) + (9«) 2 * +2 (Y') • • • (143) 

which has the form appearing in (141) as part of the hypothesis. 
Now take the case of n odd where, from (141), 



U _V 



yn +T2^2T=T-T2nF2---- (W4) 



Differentiate and rearrange to obtain, 

tf*Vi = (fa) 2i (U(2n -2i- D92) 

+ (9i) 2 ' +1 (t/' + V(2n - 2/ - 2)92) - (9i) 2 ' +2 ( V) 

As 17' < 0, V > 0, the form in (140) of the induction hypothesis is 
satisfied. This concludes the proof. 

APPENDIX C 

Proof of Proposition 10 

In the following we shall need the sign property of the derivatives of 
/(.) as given in (25) and (26). Also recall U = f'(0) = -a 
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Jo Jo 



ze~ mz) dz < ze~ Nf * z = l(Na) 2 (145) 



,-N(/ l2 +/ 2 2 2 ), 



'0 



2iv/ 2 
Similarly, 



4/2 ViV/2 

(146) 



I e~ mz) dz < -l/(aN) (147) 

Jo 

" V 4^ e02iW2(1 " er/r( " a J*W*))- < 148 > 



Thus, 



J. 


ze 
e~ 


~ Nf(z) dz 
mz) dz 


< 


y/4f 2 /ir 


e c 


2 N/4f 2 


f 


ct 2 N 3/2 [i 


- erf(- 


a VAT/4/ 2 )] 




1 + Vl + 


8f 2 /a 2 N 





Ol iw ' (149) 

where, to bound the term dependent on the error function, we have 
used the left inequality in the following 13 

(* + V* 2 + 2)" 1 < e' 2 J e~ y2 dy < [x(l + Vl + 4/(*r* 2 ))]-\ (150) 
We obtain from (146), (147) and (150) the results analogous to (149): 

ze- Nf(z) dz / \ 



f 



Nfiz) d2 



e ' az 



fe W 1 2 (151) 

2 £ \ 1 + Vl + 16f 2 /(ira 2 N) I 



Equations (149) and (151) taken together with the representation for 
the utilization given in (22) is the content of Proposition 10. 
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APPENDIX D 

Proof of Proposition 1 1 

To prove the proposition, we begin by substituting (54) in the 
expression for I in (56) to obtain after a change of variable 

7- (y/N + c) [ e- ( ^ +e " J] (1 + w/(m y/N))™**^ dp. (152) 

Jo 

We may write (152) as 

(153) 



/= f e - Ap2/2 - Bv h(u, VN)dp, 
Jo 



where 

A = YJ>i/a?, B = c- ^di/di (154) 

and 

h(v, >/N) = (>/N+c) exp[Av 2 /2 + Bv 

- p(Jn + c)][II<l + "/(a« JN)) b ' N+d -^l (155) 

The quantity h(v, y/N) has been denned so that when the second 
bracketed expression [•] is written as exp(log [•]) and then the log 
[.] term expanded, a cancellation of the leading terms is effected by 
multiplication with the exponential term in (155). In this step, notice 
has to be taken of the constraint in (55) in that exp[Av 2 /2 + Bv - 
v(y/N+ c)] = exp[Av 2 /2 - v(>/N £6,/a,- + £d,/a,-)]. In this manner we 
obtain, 

h(p, >/N) = ( JN + c)exp( £ Fj(p)/y/N J j, (156) 

where 

In (156) we have an exponential of a power series in 1/\JN which may 
equivalently be represented directly as a power series in l/\/N: 

h{v, JN) = (JN+c) £ Gj(p)/Vn j . (158) 

For example, Go = 1, Gi = F h G 2 = F 2 /2 + F 2 , G 3 = Ff/6 + F X F 2 + Fa. 
A feature to observe is that Gj(v) is a polynomial of degree 3/ in v. 
For the final form, 

h(v, s/N) = X Hjiv)/>/N J '-\ (159) 

y=o 
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where Hj(p) = Gj(p) + cGj-i(p), also a polynomial of degree 3/ in v. It 
is understood that G-\(v) = 0. 

Insertion of (159) in (153) yields Proposition 11, namely 

I=icj/JN J -\ (160) 

j-o 

where c, = f e^^^Hj^dv. (161) 

Jo 
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